rm(list = ls())
set.seed(1111)
library(MatchIt)
library(dplyr)
library(ggplot2)
library(haven)
library(stargazer)
library(tableone)
library(Hmisc)
library(cem)
library(Matching)
library(car)
library(Zelig)
library(RColorBrewer)
library(arm)
setwd("~/Dropbox/AttorneyMainData")
data <- read_dta("Attorney.dta")
data$hysLaw <- car::recode(data$SchoolLawyerArgue,"'Harvard'=1;'Yale'=1;'Stanford'=1;'Columbia'=1;'Chicago'=1;else=0")
data$demPres <- car::recode(data$SGOrals,"'McCree'=1;'Lee'=0;'(Acting) Fried'=0;'Fried'=0;'(Acting) Bryson'=0;'Starr'=0;'Days'=1;'(Acting) Dellinger'=1;'Waxman'=1;'(Acting) Underwood'=0; 'Olson'=0;'(Acting) Clement'=0;'Clement'=0;'Garre'=0;'(Acting) Kneedler'=1;'Kagan'=1;'(Acting) Katyal'=1;'Verrilli'=1;'(Acting) Gershengorn'=1;'(Acting) Francisco'=1;'(Acting) Wall'=1; 'Francisco'=1")
##Bryson was acting for both Bush 1 and Reagan
data$demPres[data$SGOrals=="(Acting) Bryson" & data$term=="1992"] <- 1
data$SGOralsCollapseF <- as.factor(data$SGOralsCollapse)
data$justiceNameF <- as.factor(data$justiceName)
data$ActingSGOrals <- as.numeric(factor(data$ActingSGOrals, ordered=F))-1
data$TopLawSchool <- as.numeric(factor(data$TopLawSchool, ordered=F))-1
data$DCFirm <- as.numeric(factor(data$DCFirm, ordered=F))-1
data$ClerkDummy <- as.numeric(factor(data$ClerkDummy, ordered=F))-1
data$OppPetOrResp <- as.numeric(factor(data$OppPetOrResp, ordered=F))-1
data$CriminalDummy <- as.numeric(factor(data$CriminalDummy, ordered=F))-1
data$CivLibDummy <- as.numeric(factor(data$CivLibDummy, ordered=F))-1
data$EconDummy <- as.numeric(factor(data$EconDummy, ordered=F))-1
data$Acela <- as.numeric(factor(data$Acela, ordered=F))-1
data$LawyerArgueOSG <- as.numeric(factor(data$LawyerArgueOSG, ordered=F))-1
data$GenderLawyerArgue <- as.numeric(factor(data$GenderLawyerArgue, ordered=F))-1
data$RaceLawyerArgue <- as.numeric(factor(data$RaceLawyerArgue, ordered=F))-1
data$IdeologyOppCounsel <- as.numeric(factor(data$IdeologyOppCounsel, ordered=F))-1
ctLevel <- subset(data, ct_level==1)
summary(data)
dim(data)


##Create Datasets
#Justice-level
xvars<-c("SGOralsCollapseF","ActingSGOrals","hysLaw","DCFirm","ClerkDummy","Rehnquist","Roberts","IdeologyOppCounsel","OppPetOrResp","demPres","justiceNameF")
data_nomiss <- data %>%  # MatchIt does not allow missing values
  dplyr::select(PriorOral, JVote, one_of(xvars)) %>%
  na.omit()
data_nomiss <- data.frame(data_nomiss)
table(data_nomiss$PriorOral)
pretab1<-CreateTableOne(vars=xvars, strata ="PriorOral", 
                        data=data_nomiss, test = FALSE)
print(pretab1, smd = TRUE)

##Case-level
xvars_case<-c("SGOralsCollapseF","ActingSGOrals","hysLaw","DCFirm","ClerkDummy","Rehnquist","Roberts","IdeologyOppCounsel","OppPetOrResp","demPres")
data_case_nomiss <- ctLevel %>%  # MatchIt does not allow missing values
  dplyr::select(PriorOral, WhoWon, minWinning, unan, nearUnan, one_of(xvars_case)) %>%
  na.omit()
data_case_nomiss <- data.frame(data_case_nomiss)
table(data_case_nomiss$PriorOral)
pretab1<-CreateTableOne(vars=xvars_case, strata ="PriorOral", 
                        data=data_case_nomiss, test = FALSE)
print(pretab1, smd = TRUE)


##CEM
match_justice_cem <- matchit(PriorOral ~ SGOralsCollapseF+ActingSGOrals+hysLaw+DCFirm+ClerkDummy+IdeologyOppCounsel+demPres+Rehnquist+Roberts+OppPetOrResp+justiceNameF,
                             method = "cem", data = data_nomiss)
summary(match_justice_cem)
dta_justice_cem <- match.data(match_justice_cem)

match_case_cem <- matchit(PriorOral ~ SGOralsCollapseF+ActingSGOrals+hysLaw+DCFirm+ClerkDummy+IdeologyOppCounsel+demPres+Rehnquist+Roberts+OppPetOrResp,
                          method = "cem", data = data_case_nomiss)
summary(match_case_cem)
dta_case_cem <- match.data(match_case_cem)

##SI, Table A3
table(data_nomiss$PriorOral)
table(data_case_nomiss$PriorOral)
table(dta_justice_cem$PriorOral)
table(dta_case_cem$PriorOral)

##Run Models
data_case_cem <- match.data(match_case_cem)
data_justice_cem <- match.data(match_justice_cem)

full_case <- glm(WhoWon~PriorOral, family="binomial", data = ctLevel); summary(full_case)
cem_case <- glm(WhoWon~PriorOral, family="binomial", data = match.data(match_case_cem), weights=match.data(match_case_cem)$weights); summary(cem_case)
full_case_m <- glm(WhoWon~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres, family="binomial", data = ctLevel); summary(full_case_m)
cem_case_m <- glm(WhoWon~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres, family="binomial", data = match.data(match_case_cem), weights=match.data(match_case_cem)$weights); summary(cem_case_m)

###SI, Table A4
stargazer(full_case,full_case_m,cem_case,cem_case_m,type="text",digits=2)

full_justice <- glm(JVote~PriorOral, family="binomial", data = data); summary(full_justice)
cem_justice <- glm(JVote~PriorOral, family="binomial", data = match.data(match_justice_cem), weights=match.data(match_justice_cem)$weights); summary(cem_justice)
full_justice_m <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+justiceNameF, family="binomial", data = data); summary(full_justice_m)
cem_justice_m <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+justiceNameF, family="binomial", data = match.data(match_justice_cem), weights=match.data(match_justice_cem)$weights); summary(cem_justice_m)

###SI, Table A5
stargazer(full_justice,full_justice_m,cem_justice,cem_justice_m,type="text",digits=2)

##simulate
set.seed(1111)
dd <- data
aa <- sim(full_justice)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,1] + acoef[,2]*1)
untreated <- invlogit(acoef[,1] + acoef[,2]*0)
j_n <- treated-untreated

dd <- data
aa <- sim(full_justice_m)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0  + acoef[,"justiceNameFJPStevens"]*1)
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0  + acoef[,"justiceNameFJPStevens"]*1)
j_m <- treated-untreated

dd <- match.data(match_justice_cem)
aa <- sim(cem_justice)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,1] + acoef[,2]*1)
untreated <- invlogit(acoef[,1] + acoef[,2]*0)
j_n_match <- treated-untreated

dd <- match.data(match_justice_cem)
aa <- sim(cem_justice_m)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0  + acoef[,"justiceNameFJPStevens"]*1)
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0  + acoef[,"justiceNameFJPStevens"]*1)
j_m_match <- treated-untreated
mean(j_m_match)

##Case-level
dd <- ctLevel; c_n <- predict(full_case,data.frame(PriorOral=c(0,1)),type="response",se=T)
dd <- ctLevel; c_m <- predict(full_case_m,data.frame(PriorOral=c(0,1),hysLaw=median(dd$hysLaw,na.omit=T),DCFirm=median(dd$DCFirm,na.omit=T),ClerkDummy=median(dd$ClerkDummy,na.omit=T),Rehnquist=median(dd$Rehnquist,na.omit=T),Roberts=median(dd$Roberts,na.omit=T),OppPetOrResp=median(dd$OppPetOrResp,na.omit=T),IdeologyOppCounsel=1,demPres=0),type="response",se=T)
dd <- match.data(match_case_cem);c_n_match <- predict(cem_case,data.frame(PriorOral=c(0,1)),type="response",se=T)
dd <- match.data(match_case_cem);c_m_match <- predict(cem_case_m,data.frame(PriorOral=c(0,1),hysLaw=median(dd$hysLaw,na.omit=T),DCFirm=median(dd$DCFirm,na.omit=T),ClerkDummy=median(dd$ClerkDummy,na.omit=T),Rehnquist=median(dd$Rehnquist,na.omit=T),Roberts=median(dd$Roberts,na.omit=T),OppPetOrResp=median(dd$OppPetOrResp,na.omit=T),IdeologyOppCounsel=1,demPres=0),type="response",se=T)

dd <- ctLevel
aa <- sim(full_case)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,1] + acoef[,2]*1)
untreated <- invlogit(acoef[,1] + acoef[,2]*0)
c_n <- treated-untreated

dd <- ctLevel
aa <- sim(full_case_m)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 )
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 )
c_m <- treated-untreated

dd <- match.data(match_case_cem)
aa <- sim(cem_case)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,1] + acoef[,2]*1)
untreated <- invlogit(acoef[,1] + acoef[,2]*0)
c_n_match <- treated-untreated

dd <- match.data(match_case_cem)
aa <- sim(cem_case_m)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 )
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0)
c_m_match <- treated-untreated

###Manuscript, Figure 5
pdf("MainResults.pdf",width=15,height=6.6,paper='special') 
par(mar=c(5.1,2, 4.1, 2.1),mfrow=c(1,2),oma=c(0,5,0,0))
plot(NA,xlim=c(-.1,.3),ylim=c(.5,2.5),ylab="",xlab="ATT", main="Justice-Level Data",yaxt="n")
axis(2, at=c(2,1),labels=c("Naive","Multivariate"), las=2)
points(y=c(2.1,1.1),x=c(mean(j_n),mean(j_m)), pch=19,cex=1.5,col="gray70")
points(y=c(1.9,.9),x=c(mean(j_n_match),mean(j_m_match)), pch=18,cex=1.75)
segments(y0=2.1,y1=2.1,x0=quantile(j_n,.025),x1=quantile(j_n,.975),lwd=2,col="gray70")
segments(y0=1.1,y1=1.1,x0=quantile(j_m,.025),x1=quantile(j_m,.975),lwd=2,col="gray70")
segments(y0=1.9,y1=1.9,x0=quantile(j_n_match,.025),x1=quantile(j_n_match,.975),lwd=2)
segments(y0=.9,y1=.9,x0=quantile(j_m_match,.025),x1=quantile(j_m_match,.975),lwd=2)
abline(v=0,lwd=3,col="gray60")

#Case-Level Figure
plot(NA,xlim=c(-.1,.3),ylim=c(.5,2.5),ylab="",xlab="ATT", main="Court-Level Data",yaxt="n")
points(y=c(2.1,1.1),x=c(mean(c_n),mean(c_m)), pch=19,cex=1.5,col="gray70")
points(y=c(1.9,.9),x=c(mean(c_n_match),mean(c_m_match)), pch=18,cex=1.75)
segments(y0=2.1,y1=2.1,x0=quantile(c_n,.025),x1=quantile(c_n,.975),lwd=2,col="gray70")
segments(y0=1.1,y1=1.1,x0=quantile(c_m,.025),x1=quantile(c_m,.975),lwd=2,col="gray70")
segments(y0=1.9,y1=1.9,x0=quantile(c_n_match,.025),x1=quantile(c_n_match,.975),lwd=2)
segments(y0=.9,y1=.9,x0=quantile(c_m_match,.025),x1=quantile(c_m_match,.975),lwd=2)
abline(v=0,lwd=3,col="gray60")
legend("bottomright",c("Full Dataset","Matched Dataset"),       col=c("gray70","black"),pch=c(19,18),pt.cex=c(1,1.3))
par(mar=c(5.1,4.1, 4.1, 2.1),oma=c(0,0,0,0))
par(mfrow=c(1,1))
dev.off()

###Manuscript, Figure 4
pdf("Coef.pdf",width=15,height=9,paper='special') 
par(mar=c(5.1,2, 4.1, 2.1),mfrow=c(1,2),oma=c(0,10,0,0))
plot(NA,xlim=c(-.75,1.25),ylim=c(.5,9.5),ylab="",xlab="Logistic Regression Coefficient", main="Justice-Level Data",yaxt="n")
axis(2, at=c(9,8,7,6,5,4,3,2,1),labels=c("Experienced Attorney","Educational Quality","Worked in DC Firm","Was a Supreme Court Clerk","Rehnquist Court","Roberts Court","Represented Petitioner","Represented Liberal Side","SG Appointed by Democrat"), las=2)
abline(v=0,lwd=3,col="gray60")
points(y=c(9,8,7,6,5,4,3,2,1)+.1,x=coef(full_justice_m)[2:10], pch=19,cex=1.5,col="gray70")
points(y=c(9,8,7,6,5,4,3,2,1)-.1,x=coef(cem_justice_m)[2:10], pch=18,cex=1.75)
for(i in 1:9){
  segments(y0=10-i+.1,y1=10-i+.1,x0=coef(full_justice_m)[i+1]+1.96*se.coef(full_justice_m)[i+1],x1=coef(full_justice_m)[i+1]-1.96*se.coef(full_justice_m)[i+1],lwd=2,col="gray70")
  segments(y0=10-i-.1,y1=10-i-.1,x0=coef(cem_justice_m)[i+1]+1.96*se.coef(cem_justice_m)[i+1],x1=coef(cem_justice_m)[i+1]-1.96*se.coef(cem_justice_m)[i+1],lwd=2)
}

plot(NA,xlim=c(-.75,1.25),ylim=c(.5,9.5),ylab="",xlab="Logistic Regression Coefficient", main="Case-Level Data",yaxt="n")
abline(v=0,lwd=3,col="gray60")
points(y=c(9,8,7,6,5,4,3,2,1)+.1,x=coef(full_case_m)[2:10], pch=19,cex=1.5,col="gray70")
points(y=c(9,8,7,6,5,4,3,2,1)-.1,x=coef(cem_case_m)[2:10], pch=18,cex=1.75)
for(i in 1:9){
  segments(y0=10-i+.1,y1=10-i+.1,x0=coef(full_case_m)[i+1]+1.96*se.coef(full_case_m)[i+1],x1=coef(full_case_m)[i+1]-1.96*se.coef(full_case_m)[i+1],lwd=2,col="gray70")
  segments(y0=10-i-.1,y1=10-i-.1,x0=coef(cem_case_m)[i+1]+1.96*se.coef(cem_case_m)[i+1],x1=coef(cem_case_m)[i+1]-1.96*se.coef(cem_case_m)[i+1],lwd=2)
}
legend("bottomright",c("Full Dataset","Matched Dataset"),       col=c("gray70","black"),pch=c(19,18),cex=.8,pt.cex=c(1,1.3))
dev.off()

###Manuscript, Figure 3
pdf("CEML1.pdf",width=15,height=6.6,paper='special') 
par(mar=c(5.1,2, 4.1, 2.1),mfrow=c(1,2),oma=c(0,10,0,0))
plot(NA,xlim=c(0,.5),ylim=c(.5,9.5),ylab="",xlab=expression(paste(L[1]," Statistic")), main="Justice-Level Data",yaxt="n")
axis(2, at=c(9,8,7,6,5,4,3,2,1),labels=c("Educational Quality","Worked in DC Firm","Was a Supreme Court Clerk","Rehnquist Court","Roberts Court","Represented Liberal Side","Represented Petitioner","SG Appointed by Democrat","Justice-Level Fixed Effects"), las=2)
points(y=c(9,8,7,6,5,4,3,2,1),x=i_justice_full$tab[c(3:11),3], pch=19,cex=1.5,col="gray70")
points(y=c(9,8,7,6,5,4,3,2,1),x=i_justice_cem$tab[c(3:11),3], pch=18,cex=1.75)

plot(NA,xlim=c(0,.5),ylim=c(.5,9.5),ylab="",xlab=expression(paste(L[1]," Statistic")), main="Court-Level Data",yaxt="n")
points(y=c(9,8,7,6,5,4,3,2,1),x=c(i_case_full$tab[c(3:10),3],100), pch=19,cex=1.5,col="gray70")
points(y=c(9,8,7,6,5,4,3,2,1),x=c(i_case_cem$tab[c(3:10),3],100), pch=18,cex=1.75)
legend("bottomright",c("Full Dataset","Matched Dataset"),       col=c("gray70","black"),pch=c(19,18),pt.cex=c(1,1.3))
dev.off()

###Manuscript, Figure 2
barData <- read_dta("CasesAllvAtty.dta")
par(mar=c(5.1,4.1, 4.1, 2.1),oma=c(0,0,0,0))
par(mfrow=c(1,1))
pdf("../CaseComparison1.pdf",width=15,height=6.6,paper='special') 
b <- barplot(as.vector(barData$FracCases),ylim=c(0,.45),xaxt="n",col=c(rep("gray30",6),rep("gray60",19),rep("gray90",13)),
             xlab="Term",ylab="Fraction of Cases Included")
b
axis(1,at=b,labels=barData$term,las=2)
legend("topright",c("Burger Court","Rehnquist Court","Roberts Court"),lty=c(1,1,1),lwd=4,col=c("gray30","gray60","gray90"),cex=.7)
dev.off()

###Manuscript, Figure 1
pointData <- read_dta("FracPriorOral.dta")
pdf("../PriorExperienceFrac1.pdf",width=9,height=6,paper='special') 
plot(x=pointData$term,y=pointData$fracprior,cex=pointData$ntotal/30,pch=19,ylim=c(.1,.75),
     xlab="Term",ylab="Fraction of Cases With Experienced Attorney")
lw1 = loess(fracprior ~ term,data=pointData)
j <- order(pointData$term)
lines(pointData$term[j],lw1$fitted[j],col="gray70",lwd=3)
#scatter.smooth(pointData$fracprior ~ pointData$term, span = 2/3, degree = 2)
dev.off()

